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A thorough understanding of the electronic structure is a necessary first step for the design of nanoelectronics, 
chemical/bio-sensors, electrocatalysts, and nanoplasmonics using graphene. As such, theoretical spectroscopic tech- 
niques to describe both direct optical excitations and collective excitations of graphene are of fundamental importance. 
Starting from density functional theory (DFT) we use the time dependent linear response within the random phase 
approximation (TDDFT-RPA) to describe the loss function — 3{e _1 (q, us)} for graphene. To ensure any spurious in- 
teractions between layers are neglected, we employ both a radial cutoff of the Coulomb kernel, and extra vacuum 
directly at the TDDFT-RPA level. 



1 Introduction Understanding of the momentum de- 
pendent response function is one of the major challenges in 
solid state spectroscopy. One typical example is the analy- 
sis of the two particle excitation spectra as determined by 
the loss function of inelastic electron scattering. In con- 
trast to direct optical excitations, such electron energy loss 
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Figure 1 Schematic of the or- 
thorhombic graphene unit cell 
repeated twice in the surface plane. 
The x-direction corresponds to the 
zigzag direction or circumference of a 
zigzag SWNT, while the y-direction 
corresponds to the armchair direction 
or circumference of an armchair 
SWNT. The ^-direction is normal to 
the graphene surface. 



spectroscopy gives a direct probe of both collective longi- 
tudinal excitations (inter and intra-band plasmons), as well 
as single particle excitations. The momentum dependence 
of these plasmon excitations as well as single particle ex- 
citation gives a measure of the combined dispersion of the 
valence and conduction bands. Additionally, one can ex- 
tract information on the screening and dispersion of free 
charge carriers in low dimensional systems. 

To address these questions, we have performed den- 
sity functional theory (DFT) calculations within the pro- 
jector augmented wavefunction (PAW) method, and em- 
ployed the Casida method to calculate the non-interacting 
density-density response function within the random-phase 
approximation (TDDFT-RPA), to obtain the loss function 
— S{e~ 1 (q,w) [1,2]. Using this methodology, we have 
calculated the loss function for graphene along the arm- 
chair or y-direction, as shown schematically in Fig. 1 . 

2 Methodology All DFT calculations were per- 
formed using the real-space projector augmented wave- 
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function (PAW) method code GPAW [3], with a grid spacing 
of 0.2 A, and the PBE exchange correlation (xc)-functional 
[4]. An electronic temperature of kgT » 0.05 eV was used 
to obtain the occupation of the Kohn-Sham orbitals, with 
all energies extrapolated to T — K, and one unoccupied 
band per C atom included to improve convergence. 

Structural minimization was performed within the 
Atomic Simulation Environment (ASE)[5], until a maxi- 
mum force below 0.05 eV/A was obtained. An orthorhom- 
bic 2.46 A x 4.26 A x 8.00 A supercell was employed, 
consisting of four C atoms, as depicted in Fig. 1. Non- 
periodic boundary conditions were enforced in the in- 
direction normal to the graphene surface, so that both 
the electron density and Kohn-Sham wavefunctions — > 
as z — > or z — > L, where L is the length of the unit cell 
in the z -direction. A Monkhorst-Pack fc-point sampling of 
25 fc-points along the zigzag direction, and 15 fc-points 
along the armchair direction of the graphene surface was 
employed, yielding a longitudinal momentum transfer res- 
olution Aq of 0.102 A -1 and 0.098 A -1 respectively. 
14 unoccupied bands per C atom were converged in the 
self-consistent calculation, which was ultimately found to 
be more than sufficient to converge the loss function for 
energies up to 50 eV, as shown in Fig. 2. 

Calculations of the loss function have been performed 
within the Casida methodology [1], employing time de- 
pendent density functional theory within the random 
phase approximation (TDDFT-RPA), as recently imple- 
mented within GPAW [6]. Within this framework the non- 
interacting density-density function Xqg' (l> w ) ^ or mo ~ 
mentum transfer q at energy u is given by 
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Here the sum is over reciprocal lattice vectors k and band 
numbers n and n', with e„k the eigenenergy of the n th band 
at k, /„k the Fermi-Dirac occupation of the n th band at 
k, 7 the peak broadening, fl the volume of the supercell, 
G and G' the reciprocal unit cell vectors, and ^> n k(r) the 
real-space Kohn-Sham wavefunctions for the n th band with 
reciprocal lattice-vector k. 

Including local field effects, we may write the in- 
verse macroscopic dielectric function £ _1 (q, u)) within the 
random phase approximation (RPA) in terms of the non- 
interacting density-density response function Xgg'CI'^) 
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where 5qq' is the Kronecker delta, and vq (q) is the Fourier 
transform of the Coulomb kernel. It should be noted that 
the inclusion of exchange and correlation effects in vq (q) 



at the LDA level adds a minor correction to the present 
results, as already shown for the case of graphite [7,8]. 

As discussed in Ref. [9], for a 3D periodic system with 
translational invariance, the Coulomb kernel is 
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However, for a system which is periodic in only two di- 
mensions, such as a bulk slab or graphene, interactions be- 
tween periodic images in a TDDFT-RPA calculation may 
be significant due to the long-range behaviour of v 3D . This 
will be the case even for systems with sufficient vacuum 
to converge the electron density at the DFT level. On the 
other hand, image — image interactions are included at the 
TDDFT-RPA level only through v 3D . This motivates us to 
introduce a 2D periodic Coulomb kernel, v 2D , which is 
both translationally invariant and zero for \z\ > R, where 
R is the "radial cutoff" for the Coulomb kernel. In this case 
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Employing the suggested choice of R — L/2 [9], since 
Gz =n|-, where n G Z, we find 
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From Eqn. (4) we clearly see that for L S> 2/q or q > 
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Alternatively, we may introduce further regions of vac- 
uum separating the images directly at the TDDFT-RPA 
level. Since in the added vacuum regions the matrix el- 
ements for the occupied Kohn-Sham wavefunctions are 
zero, the inclusion of extra vacuum in Eqn. (1) only enters 
into the non-interacting density-density response function 
through the unit cell volume Q and the reciprocal unit cell 
vectors G. We may thus introduce extra unit cells of vac- 
uum, or "zero padding" in the non-periodic direction, by 
doubling or tripling L when computing the set of G vectors 
to include at the TDDFT-RPA level. In this way, increas- 
ing the length of the unit cell in the non-periodic direction 
through the inclusion of vacuum effectively increases the 
density of sampling of the reciprocal unit cell. 

Finally, the quantities of fundamental interest are the 
loss function — 3{e _1 (q, uj)}, the adsorption or imaginary 
part of the dielectric function 3{e(q, w)}, and the real part 
of the dielectric function 3?{e(q,w)}, which may be ob- 
tained from Eqn. (2). 

3 Results & Discussion To test the convergence of 
the TDDFT-RPA calculated dielectric function with re- 
spect to the number of unoccupied bands included within 
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Figure 2 Convergence with respect to the number of 
unoccupied bands n unocc of the calculated loss function 
— 3{e~ 1 (q,a;)}, imaginary part of the dielectric func- 
tion 3{e(q,u;)}, and real part of the dielectric function 
5ft{e(q, u>)} versus energy uj in eV for momentum transfer 
||q|| fsO.l A -1 along the armchair direction of graphene. 
Contributions from including ri unocc = 1, 2, 3, 4, 5, 6, 7, 
8, 10, and 14 unoccupied bands per C atom are shown. A 
plane-wave cutoff energy e cut = 60 eV, corresponding to 89 
G vectors was used. 

Table 1 Percentage error for the /-sum rule Af, and con- 
vergence with the number of unoccupied bands n unocc per 
C atom. 
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the calculation n unocc , we make use of the /-sum rule, 

j ^9{e( w )} = -^, (5) 
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Figure 3 Influence of neighbouring image interactions 
on the calculated loss function — 3{e _1 (q, uj)}, imaginary 
part of the dielectric function 3{e(q, u>)}, and real part of 
the dielectric function 5R{e(q,a;)} versus energy uj in eV 
for momentum transfer ||q|| « 0.1 A -1 along the armchair 
direction of graphene, and including n unocc = 8 bands/atom. 
Results are shown for e cut = 15 eV with 8 A of vacuum ( — 

— ), 16 A of vacuum (- • -), and 24 A of vacuum ( ), 

corresponding to 1 1, 25, and 37 G vectors respectively, and 
for e cut = 60 eV with 8 A of vacuum, corresponding to 89 

G vectors, both with ( ) and without ( ) applying a 

radial cutoff R = 4 A to the Coulomb kernel. 



where N e is the number of electrons in the unit cell of vol- 
ume SI. The percentage error in the adherence of the imag- 
inary part of the calculated dielectric function to the sum 
rule, Af, for n unocc = 1 to 14 unoccupied bands per C 
atom is shown in Table 1 . For n unocc = 8 bands/atom, the 
/-sum rule is quite well satisfied, considering that the inte- 
gral has been truncated at cj max = 60 eV, and the approxi- 
mation du> w Auj = 0.02 eV. It should also be noted that 
the /-sum rule is not satisfied exactly when using non-local 
pseudopotentials, although the correction is within the er- 
rors given in Table 1 [2] . 

As a more direct test, we plot in Fig. 2 the loss func- 
tion and both the real and imaginary parts of the dielectric 
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Figure 4 Convergence with respect to the plane- 
wave cutoff energy e cut for the calculated loss func- 
tion — Q ! {£~ 1 (q, oj)} versus energy u> in eV and momen- 
tum transfer ||q|| in A -1 along the armchair direction of 
graphene, including n unocc = 8 bands/atom, and a radial 

cutoff R = 4 A. Results for e cut = 15 ( ), 30 (- 

-), 60 eV ( ), and 90 eV ( ), corresponding to 11, 

27, 89, and 149 G vectors respectively, are shown. 



function for n unoC c between 1 and 14 bands per C atom. We 
consider the limit of low momentum (||q|| « 0. 1 A~ 1 in the 
armchair direction) since in this limit Xqc w ^ ^ e rat her 
insensitive to number of G vectors included. 

Figure 2 clearly shows that (1) the dielectric func- 
tion converges faster than the loss function with respect 
to n unocc , (2) for n unocc = 8 bands/atom we have conver- 
gence up to about 50 eV, and (3) with only 2-3 unoccu- 
pied bands per C atom the loss function is converged semi- 
quantitatively up to 10 eV, as shown in the inset. Combined 
with the values for the /-sum rule error in Table 1, this sug- 
gests that n unocc = 8 bands/atom is sufficient to describe 
the dielectric response of graphene-like materials up to 50 
eV. Further, a near quantitative description of the behaviour 
of the low energy ir plasmons is obtained by including only 
2-3 unoccupied bands per C atom. 

However, although both the electron density and Kohn- 
Sham wavefunction have been set to zero at z = and 
z = L through the use of non-periodic boundary condi- 
tions at the DFT level, interactions between periodic im- 
ages in the z-direction are included through the Coulomb 
kernel in Eqn. (2). Figure 3 shows how much neighbour- 
ing image interactions contribute to the loss and dielectric 
functions for low momentum transfers. By either applying 
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Figure 5 Intensity of the calculated loss function 
— S{e~ 1 (q,w)} as a function of energy u in eV and mo- 
mentum transfer ||q|| in A -1 along the armchair direction 
of graphene, with n unocc = 8 bands/atom, e cut = 60 eV, cor- 
responding to 89 G vectors, and using a radial cutoff R = 
4 A. 



a radial cutoff of R = 4 A using v 2D or increasing the unit 
cell length L to 24 A at the TDDFT-RPA level, we ob- 
tain the same converged loss and dielectric functions, with 
image — image interactions removed. However, as the addi- 
tion of extra unit cells of vacuum, even at only the TDDFT- 
RPA level, increases the size of the Xqq' matrix, we shall 
employ a radial cutoff of R = 4 A from hereon. 

Figure 4 shows the convergence of the loss function 
dispersion with respect to the plane-wave energy cutoff 
e cu t, and the corresponding number of G vectors included 
in the calculation. As mentioned above, for low to medium 
momentum transfers (||q|| < 0.8 A -1 ) only a few G vec- 
tors (10 — 30) are required to converge the loss function. 
Further, even at high momentum transfers (||q|| > 1 A -1 ), 
the differences between e cut = 15 and 90 eV lie mainly in 
the intensities rather than the peak positions. This is espe- 
cially true if we restrict consideration to the perhaps most 
interesting low energy ir plasmon regime. 

The intensity of the converged loss function for 
graphene is plotted in Fig. 5 for n unocc = 8 bands/atom, 89 
G vectors, and a radial cutoff R = 4 A. We clearly see the 
near linearly dispersive ir plasmon at 5 — 10 eV, the non- 
dispersive ir plasmon at about 5.5 eV, and the dispersing 
a + ir plasmon between 15-30 eV, from experiment [8]. 
We also find that the optical limit ||q|| — > + for graphene 
is well reproduced [10]. 
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4 Conclusions In conclusion, we have extended the 
TDDFT-RPA implementation within GPAW to employ both 
a radial cutoff of the Coulomb kernel v 2D for 2D peri- 
odic systems, and include extra unit cells of vacuum at 
the TDDFT-RPA level. We find these spurious image — 
image interactions have a significant impact on the calcu- 
lated loss function for isolated systems at low momentum 
transfer (||q|| < 0.5 A -1 ), as demonstrated for graphene. 
We also find for carbon systems, including about 3 unoc- 
cupied bands per atom is sufficient to obtain a near quanti- 
tative description of the loss function up to 10 eV. 

These results are particularly important in the area of 
nanoplasmonics, for the description of the low energy free- 
charge carrier plasmons induced by electrostatic or potas- 
sium doping. Such systems can be well described within a 
simple RPA model by shifting the Fermi level rigidly, as 
we will show in future work. 
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